options(scipen = 999, digits = 4)
knitr::opts_chunk$set(comment = "#")
r = getOption("repos")
r["CRAN"] = "http://cran.us.r-project.org"

options(repos = r)

## load required packages
ipak <- function (pkg) {
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}

packages <- c("tidyverse", "Hmisc", "lattice", "multcomp", "lsmeans", "schoRsch", "influence.ME", "lme4", "effects", "lmerTest", "cowplot", "irr", "simr", "plyr", "dplyr","patchwork", "wesanderson", "MuMIn", "devtools", "dplyr", "ggResidpanel", "HLMdiag", "mixed", "sjPlot")

ipak(packages)

## Warning: package 'mixed' is not available (for R version 4.0.2)

## Warning: replacing previous import 'vctrs::data_frame' by 'tibble::data_frame'
## when loading 'dplyr'

## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'mixed'

##    tidyverse        Hmisc      lattice     multcomp      lsmeans     schoRsch 
##         TRUE         TRUE         TRUE         TRUE         TRUE         TRUE 
## influence.ME         lme4      effects     lmerTest      cowplot          irr 
##         TRUE         TRUE         TRUE         TRUE         TRUE         TRUE 
##         simr         plyr        dplyr    patchwork  wesanderson        MuMIn 
##         TRUE         TRUE         TRUE         TRUE         TRUE         TRUE 
##     devtools        dplyr ggResidpanel      HLMdiag        mixed       sjPlot 
##         TRUE         TRUE         TRUE         TRUE        FALSE         TRUE

# set summed contrasts
options(contrasts = c("contr.sum", "contr.poly")) 

# fix bs with dplyr 
detach("package:dplyr", unload = TRUE)

## Warning: 'dplyr' namespace cannot be unloaded:
##   namespace 'dplyr' is imported by 'tidyr', 'sjPlot', 'plotly', 'broom', 'janitor', 'dbplyr', 'sjmisc', 'sjstats', 'qqplotr', 'HLMdiag' so cannot be unloaded

library("dplyr")

ind.data.pre <- read.csv("./individual_data.csv") %>%
  mutate(experiment=str_sub(condition,1,1)) %>%
  mutate_if(is.character,as.factor) %>%
  mutate(subj = as.factor(subj))%>%
  mutate_at(.vars="condition", .funs = tolower)

ma.data <- read.csv("./ma_data.csv") %>%
  rename(paper = study_ID,
         condition = expt_condition,
         experiment = expt_num,
         task = looking_task) %>%
    mutate_if(is.character,as.factor) %>%
  mutate_at(.vars="condition", .funs = tolower)

study.design <- ma.data %>%
  select(paper, experiment, condition, task, training_yesno, action_consequence, actor_hand, agent_efficient_fam, object_diff_size_huge, action_causal, action_consequence, location_object_goal_ambiguous, bothobjects_present_visible_fam, agent, PI_group) %>%
    mutate_if(is.character,as.factor) %>%
  mutate(experiment = as.factor(experiment)) %>%
  filter(paper != "sommerville2005",
         task != "causes")


# merge study design with individual looks from babies

ind.data <- left_join(ind.data.pre, study.design, 
                      by=c("paper", "experiment", "condition")) %>%
  mutate(sex = tolower(str_sub(sex,1,1)),
         look_pref = unexp_look - exp_look) %>%
  mutate(paper = str_replace_all(paper, "unpublished", "unpub")) %>%
  mutate(task = str_replace_all(task, "efficiency", "constraints")) %>%
  mutate_if(is.character,as.factor) 

str(ind.data)

# 'data.frame': 626 obs. of  20 variables:
#  $ subj                           : Factor w/ 626 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
#  $ condition                      : Factor w/ 25 levels "1_active","1_control",..: 12 12 12 12 12 12 12 12 12 12 ...
#  $ ageday                         : num  103.3 86.9 105.3 129.8 107.3 ...
#  $ exp_look                       : num  13.1 18.4 15.1 7.1 5.8 9 18.9 6 31 42.9 ...
#  $ unexp_look                     : num  60 60 46.3 13.6 12 60 8 21 60 42.7 ...
#  $ sex                            : Factor w/ 2 levels "f","m": 2 2 2 1 1 1 2 2 2 1 ...
#  $ paper                          : Factor w/ 8 levels "choi_unpub","choi2018",..: 6 6 6 6 6 6 6 6 6 6 ...
#  $ experiment                     : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
#  $ task                           : Factor w/ 2 levels "constraints",..: 2 2 2 2 2 2 2 2 2 2 ...
#  $ training_yesno                 : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
#  $ action_consequence             : Factor w/ 3 levels "location_change",..: 2 2 2 2 2 2 2 2 2 2 ...
#  $ actor_hand                     : Factor w/ 3 levels "bare","gloved",..: NA NA NA NA NA NA NA NA NA NA ...
#  $ agent_efficient_fam            : Factor w/ 2 levels "no","yes": NA NA NA NA NA NA NA NA NA NA ...
#  $ object_diff_size_huge          : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
#  $ action_causal                  : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
#  $ location_object_goal_ambiguous : Factor w/ 2 levels "no","yes": 2 2 2 2 2 2 2 2 2 2 ...
#  $ bothobjects_present_visible_fam: Factor w/ 3 levels "","no","yes": 3 3 3 3 3 3 3 3 3 3 ...
#  $ agent                          : Factor w/ 3 levels "hand","object",..: 2 2 2 2 2 2 2 2 2 2 ...
#  $ PI_group                       : Factor w/ 5 levels "desrochers","luo",..: 2 2 2 2 2 2 2 2 2 2 ...
#  $ look_pref                      : num  46.9 41.6 31.2 6.5 6.2 ...

# View(ind.data)

# function that returns column of standardized betas from lmer model
gen.beta <- function(model) {
    df <- data.frame(fixef(model))
    names(df) <- c("beta")
    return(df)
}

# function that computes CIs and returns them in df
gen.ci <- function(model) {
  df <- data.frame(confint(model))
  names(df) <- c("lower", "upper")
  return(df)
}

# function that converts model summary (lmer) to df
gen.m <- function(model) {
  df <- data.frame(coef(summary(model)))
  names(df) <- c("b", "se", "df", "t", "p")
  return(df)
}

# function that converts model summary (lm) to df
gen.lm <- function(model) {
  df <- data.frame(coef(summary(model)))
  names(df) <- c("b", "se", "t", "p")
  return(df)
}

# function that returns age info and number of females in a dataset
ages <- function(longdata) {
  longdata %>% summarize(mean = mean(ageday), min=range(ageday)[1], max=range(ageday)[2], f=sum(sex=="f")/2)
}

# function that returns formatted result from lme4/lmerTest table
report <- function(table, index, places, tails, flip) {
  if (tails == "1") {
    p <- round(table$p[index], places)/2
    howmanytails <- "one-tailed"
  } else {
    p <- round(table$p[index], places)
    howmanytails <- "two-tailed"
  }
  if (p < .001) {
    p <- "<.001"
  } else {
    p <- paste("=", round(p, places), sep = "")
  }
  if (missing(flip)) {
    result <- paste("[", round(table$lower[index], places), ",", round(table$upper[index], places), "], ß=", round(table$beta[index], places), ", B=", round(table$b[index],places), ", SE=", round(table$se[index],places), ", p", p, ", ", howmanytails, sep = "")
  } else {
    result <- paste("[", -round(table$upper[index], places), ",", -round(table$lower[index], places), "], ß=", -round(table$beta[index], places), ", B=", -round(table$b[index],places), ", SE=", round(table$se[index],places), ", p", p, ", ", howmanytails, sep = "")
  }
  return(result)
}

describe <- function(dataset){
  summary <- dataset %>%
  summarise(unexp_avg = mean(unexp_look),
            exp_avg = mean(exp_look),
            lookdiff_avg = mean(look_pref),
            n = n())
  paste(
        #"M_unexp = ", summary$unexp_avg, "s ",
        #"M_exp = ", summary$exp_avg, "s ",
        "Mean_diff = ", round(summary$lookdiff_avg,3), "seconds"
  )
}

## Retrieved from : http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/#error-bars-for-within-subjects-variables
## Gives count, mean, standard deviation, standard error of the mean, and confidence interval (default 95%).
##   data: a data frame.
##   measurevar: the name of a column that contains the variable to be summariezed
##   groupvars: a vector containing names of columns that contain grouping variables
##   na.rm: a boolean that indicates whether to ignore NA's
##   conf.interval: the percent range of the confidence interval (default is 95%)
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
  library(plyr)
  
  # New version of length which can handle NA's: if na.rm==T, don't count them
  length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)
  }
  
  # This does the summary. For each group's data frame, return a vector with
  # N, mean, and sd
  datac <- ddply(data, groupvars, .drop=.drop,
                 .fun = function(xx, col) {
                   c(N    = length2(xx[[col]], na.rm=na.rm),
                     mean = mean   (xx[[col]], na.rm=na.rm),
                     sd   = sd     (xx[[col]], na.rm=na.rm)
                   )
                 },
                 measurevar
  )
  
  # Rename the "mean" column    
  datac <- plyr::rename(datac, c("mean" = measurevar))
  
  datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
  
  # Confidence interval multiplier for standard error
  # Calculate t-statistic for confidence interval: 
  # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  
  return(datac)
}
## Norms the data within specified groups in a data frame; it normalizes each
## subject (identified by idvar) so that they have the same mean, within each group
## specified by betweenvars.
##   data: a data frame.
##   idvar: the name of a column that identifies each subject (or matched subjects)
##   measurevar: the name of a column that contains the variable to be summariezed
##   betweenvars: a vector containing names of columns that are between-subjects variables
##   na.rm: a boolean that indicates whether to ignore NA's
normDataWithin <- function(data=NULL, idvar, measurevar, betweenvars=NULL,
                           na.rm=FALSE, .drop=TRUE) {
  library(plyr)
  
  # Measure var on left, idvar + between vars on right of formula.
  data.subjMean <- ddply(data, c(idvar, betweenvars), .drop=.drop,
                         .fun = function(xx, col, na.rm) {
                           c(subjMean = mean(xx[,col], na.rm=na.rm))
                         },
                         measurevar,
                         na.rm
  )
  
  # Put the subject means with original data
  data <- merge(data, data.subjMean)
  
  # Get the normalized data in a new column
  measureNormedVar <- paste(measurevar, "_norm", sep="")
  data[,measureNormedVar] <- data[,measurevar] - data[,"subjMean"] +
    mean(data[,measurevar], na.rm=na.rm)
  
  # Remove this subject mean column
  data$subjMean <- NULL
  
  return(data)
}

## Summarizes data, handling within-subjects variables by removing inter-subject variability.
## It will still work if there are no within-S variables.
## Gives count, un-normed mean, normed mean (with same between-group mean),
##   standard deviation, standard error of the mean, and confidence interval.
## If there are within-subject variables, calculate adjusted values using method from Morey (2008).
##   data: a data frame.
##   measurevar: the name of a column that contains the variable to be summariezed
##   betweenvars: a vector containing names of columns that are between-subjects variables
##   withinvars: a vector containing names of columns that are within-subjects variables
##   idvar: the name of a column that identifies each subject (or matched subjects)
##   na.rm: a boolean that indicates whether to ignore NA's
##   conf.interval: the percent range of the confidence interval (default is 95%)
summarySEwithin <- function(data=NULL, measurevar, betweenvars=NULL, withinvars=NULL,
                            idvar=NULL, na.rm=FALSE, conf.interval=.95, .drop=TRUE) {
  
  # Ensure that the betweenvars and withinvars are factors
  factorvars <- vapply(data[, c(betweenvars, withinvars), drop=FALSE],
                       FUN=is.factor, FUN.VALUE=logical(1))
  
  if (!all(factorvars)) {
    nonfactorvars <- names(factorvars)[!factorvars]
    message("Automatically converting the following non-factors to factors: ",
            paste(nonfactorvars, collapse = ", "))
    data[nonfactorvars] <- lapply(data[nonfactorvars], factor)
  }
  
  # Get the means from the un-normed data
  datac <- summarySE(data, measurevar, groupvars=c(betweenvars, withinvars),
                     na.rm=na.rm, conf.interval=conf.interval, .drop=.drop)
  
  # Drop all the unused columns (these will be calculated with normed data)
  datac$sd <- NULL
  datac$se <- NULL
  datac$ci <- NULL
  
  # Norm each subject's data
  ndata <- normDataWithin(data, idvar, measurevar, betweenvars, na.rm, .drop=.drop)
  
  # This is the name of the new column
  measurevar_n <- paste(measurevar, "_norm", sep="")
  
  # Collapse the normed data - now we can treat between and within vars the same
  ndatac <- summarySE(ndata, measurevar_n, groupvars=c(betweenvars, withinvars),
                      na.rm=na.rm, conf.interval=conf.interval, .drop=.drop)
  
  # Apply correction from Morey (2008) to the standard error and confidence interval
  #  Get the product of the number of conditions of within-S variables
  nWithinGroups    <- prod(vapply(ndatac[,withinvars, drop=FALSE], FUN=nlevels,
                                  FUN.VALUE=numeric(1)))
  correctionFactor <- sqrt( nWithinGroups / (nWithinGroups-1) )
  
  # Apply the correction factor
  ndatac$sd <- ndatac$sd * correctionFactor
  ndatac$se <- ndatac$se * correctionFactor
  ndatac$ci <- ndatac$ci * correctionFactor
  
  # Combine the un-normed means with the normed results
  merge(datac, ndatac)
}

ind.data.long <- ind.data %>%
         rename(expected = exp_look,
                unexpected = unexp_look) %>%
         pivot_longer(cols = c(expected, unexpected),
                      names_to = "trial_type",
                      values_to = "looking_time")

ind.data.summary <- summarySEwithin(data = ind.data.long, measurevar = "looking_time", betweenvars = c("task", "paper", "experiment", "condition"), withinvars = "trial_type")

# Automatically converting the following non-factors to factors: trial_type

n_infants <- ind.data.summary %>%
  filter(trial_type == "expected") %>% # just take one row per study
  summarise(sum(N))

n_conditions <- ind.data.summary %>%
  filter(trial_type == "expected") %>% # just take one row per study
  nrow()

n_papers <- ind.data.summary %>%
  filter(trial_type == "expected") %>% # just take one row per study
  group_by(paper) %>%
  count(paper) %>%
  nrow()

Summary of data sources

We searched for all journal papers, theses, and conference proceedings that reported looking time data from typically developing infants between 2 and 4 months of age in a task structured like the goals or constraints task. We also emailed two listservs (Cognitive Development Society, and Infant Studies) to request more datasets. This resulted in a final list of 11 papers and 35 conditions. We contacted the authors and asked them to send us the original datasets from this past published and unpublished work. We were able to gather original datasets from 8 papers, 29 conditions, and 626 infants (age range 72-137. A team of researchers were then randomly assigned to look through relevant papers including supplemental materials (with SL double-checking every entry and resolving disagreements). When exact values were not provided, or when we came across other ambiguities (e.g. numbers reported in the paper differ from the numbers calculated using the raw data), the team contacted authors to try and address, and also used the tool WebplotDigitizer (https://automeris.io/WebPlotDigitizer/) to extract estimated values from figures. If there were discrepancies between the paper, figures and/or raw data, we prioritized the values from the raw data if available, then author correspondence, then paper, then estimates from figures, in that order. For individual datasets, authors were asked to provide their data and a codebook, and were asked for permission to share their stimuli and data publicly on OSF. Of 8 papers (29, conditions), data from 25 conditions and stimuli from 13 conditions (either actual study videos, or example stimuli) are publicly available at https://osf.io/zwncg/.

table.goals <- ind.data %>%
  filter(task=="goals") %>%
  group_by(paper, condition, training_yesno, action_causal, action_consequence, location_object_goal_ambiguous, bothobjects_present_visible_fam, agent) %>%
  summarise(n = n(), agemin = floor(range(ageday)[1]), agemax = floor(range(ageday)[2])) %>%
  select(paper, condition, n, agemin, agemax, training_yesno, action_causal, action_consequence, location_object_goal_ambiguous, bothobjects_present_visible_fam, agent)

# `summarise()` regrouping output by 'paper', 'condition', 'training_yesno', 'action_causal', 'action_consequence', 'location_object_goal_ambiguous', 'bothobjects_present_visible_fam' (override with `.groups` argument)

table.constraints <- ind.data %>%
  filter(task=="constraints") %>%
  group_by(paper, condition, training_yesno, action_causal, action_consequence, location_object_goal_ambiguous, agent_efficient_fam, actor_hand) %>%
  summarise(n = n(), agemin = floor(range(ageday)[1]), agemax = floor(range(ageday)[2])) %>%
  select(paper, condition, n, agemin, agemax, training_yesno, action_causal, action_consequence, agent_efficient_fam, actor_hand)

# `summarise()` regrouping output by 'paper', 'condition', 'training_yesno', 'action_causal', 'action_consequence', 'location_object_goal_ambiguous', 'agent_efficient_fam' (override with `.groups` argument)

# Adding missing grouping variables: `location_object_goal_ambiguous`

knitr::kable(table.goals)

paper	condition	n	agemin	agemax	training_yesno	action_causal	action_consequence	location_object_goal_ambiguous	bothobjects_present_visible_fam	agent
choi_unpub	1_equidistant	24	76	137	no	no	none	yes	yes	person
choi_unpub	1_fartarget	24	72	130	no	no	none	yes	no	person
choi_unpub	1_neartarget	24	74	127	no	no	none	yes	yes	person
choi2018	1_hidden	16	75	121	no	no	none	yes	yes	person
choi2018	1_oneobject	16	75	127	no	no	none	yes	yes	person
choi2018	1_twoobject	16	77	130	no	no	none	yes	yes	person
gerson2014a	1_active	24	91	121	yes	no	none	yes	yes	hand
gerson2014a	1_control	24	91	120	no	no	none	yes	yes	hand
gerson2014a	1_observation	24	91	121	no	no	none	yes	yes	hand
gerson2014b	1_active	30	91	125	yes	no	none	yes	yes	hand
gerson2014b	1_observation	30	92	125	no	no	none	yes	yes	hand
gerson2014b	2_generalization	30	94	125	yes	no	none	yes	yes	hand
luo2011	1_oneobject	12	76	124	no	no	none	yes	no	object
luo2011	1_twoobject	12	79	129	no	no	none	yes	yes	object
luo2011	2_differentpositions	12	81	118	no	no	none	no	no	object
woo_unpub	1_statechange_objectswap	20	93	120	no	yes	state_change	yes	yes	person
woo_unpub	2_disambiguatingobjectgoal	24	91	121	no	yes	state_change	no	yes	person

knitr::kable(table.constraints)

location_object_goal_ambiguous	paper	condition	n	agemin	agemax	training_yesno	action_causal	action_consequence	agent_efficient_fam	actor_hand
NA	liu2019	1_pickupglove	20	92	122	no	yes	location_change	yes	gloved
NA	liu2019	2_pickupbarehand	20	93	120	no	yes	location_change	yes	bare
NA	liu2019	3_statechangebarrier	20	91	122	no	yes	state_change	yes	gloved
NA	liu2019	3_statechangenobarrier	20	91	122	no	yes	state_change	no	gloved
NA	liu2019	4_statechangenotcausal	20	93	121	no	no	state_change	yes	gloved
NA	liu2019	5_statechangecausal	26	92	121	no	yes	state_change	yes	gloved
NA	liu2019	5_statechangenotcausal	26	93	120	no	no	state_change	yes	gloved
NA	skerry2013	1_effectiveactiontraining	20	93	121	yes	yes	location_change	yes	mittened
NA	skerry2013	2_ineffectiveactiontraining	20	93	122	no	yes	location_change	yes	mittened
NA	skerry2013	3_notraining	20	91	121	no	yes	location_change	yes	mittened
NA	skerry2013	4_constrainedactionhabituation	26	93	120	yes	yes	location_change	yes	mittened
NA	skerry2013	5_unconstrainedactionhabituation	26	93	122	yes	yes	location_change	no	mittened

str(ind.data)

# 'data.frame': 626 obs. of  20 variables:
#  $ subj                           : Factor w/ 626 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
#  $ condition                      : Factor w/ 25 levels "1_active","1_control",..: 12 12 12 12 12 12 12 12 12 12 ...
#  $ ageday                         : num  103.3 86.9 105.3 129.8 107.3 ...
#  $ exp_look                       : num  13.1 18.4 15.1 7.1 5.8 9 18.9 6 31 42.9 ...
#  $ unexp_look                     : num  60 60 46.3 13.6 12 60 8 21 60 42.7 ...
#  $ sex                            : Factor w/ 2 levels "f","m": 2 2 2 1 1 1 2 2 2 1 ...
#  $ paper                          : Factor w/ 8 levels "choi_unpub","choi2018",..: 6 6 6 6 6 6 6 6 6 6 ...
#  $ experiment                     : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
#  $ task                           : Factor w/ 2 levels "constraints",..: 2 2 2 2 2 2 2 2 2 2 ...
#  $ training_yesno                 : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
#  $ action_consequence             : Factor w/ 3 levels "location_change",..: 2 2 2 2 2 2 2 2 2 2 ...
#  $ actor_hand                     : Factor w/ 3 levels "bare","gloved",..: NA NA NA NA NA NA NA NA NA NA ...
#  $ agent_efficient_fam            : Factor w/ 2 levels "no","yes": NA NA NA NA NA NA NA NA NA NA ...
#  $ object_diff_size_huge          : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
#  $ action_causal                  : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
#  $ location_object_goal_ambiguous : Factor w/ 2 levels "no","yes": 2 2 2 2 2 2 2 2 2 2 ...
#  $ bothobjects_present_visible_fam: Factor w/ 3 levels "","no","yes": 3 3 3 3 3 3 3 3 3 3 ...
#  $ agent                          : Factor w/ 3 levels "hand","object",..: 2 2 2 2 2 2 2 2 2 2 ...
#  $ PI_group                       : Factor w/ 5 levels "desrochers","luo",..: 2 2 2 2 2 2 2 2 2 2 ...
#  $ look_pref                      : num  46.9 41.6 31.2 6.5 6.2 ...

ggplot(ind.data, aes(x=look_pref)) +
  geom_histogram()+
  facet_grid(~paper)+
  theme_cowplot(12)

# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Plots

plot1 <- ggplot(ind.data.long,
       aes(trial_type, looking_time, colour=paper))+
  theme_cowplot(12)+
  theme(legend.title = element_blank())+
  geom_boxplot(outlier.alpha = 0.2)+
  geom_point(alpha=0.2)+
  geom_line(aes(group=subj), alpha=0.2)+
  stat_summary(geom="point", fun="mean", colour="black")+
  # geom_errorbar(data = ind.data.summary, colour="red", position = position_dodge(width = 5), width = 0, aes(ymin=looking_time-se, ymax=looking_time+se)) +
  stat_summary(geom="errorbar", fun.data="mean_se", width=0.2, colour="black")+
  ylab("Looking Time (s)")+
  xlab("Condition")+
  facet_wrap(~task+paper+condition, nrow=4, labeller=label_wrap_gen())+
  theme(axis.text.x = element_text(angle = 90, hjust=0.95, vjust =0.2),
        legend.position="bottom")
  
plot1

plot2.goals <- ggplot(ind.data %>% filter(task == "goals"),
                 aes(condition, look_pref, fill=paper)) +
  geom_boxplot(outlier.alpha = 0.2, position = position_dodge2(preserve = "single"))+
  geom_hline(yintercept = 0)+
  geom_point(alpha=0.2)+
  stat_summary(geom="point", fun="mean", colour="black")+
  stat_summary(geom="errorbar", fun.data="mean_se", width=0.2, colour="black")+
  theme_cowplot(18)+
  facet_grid(~paper, scales = "free_x",space = "free", labeller = label_wrap_gen())+
  ylab("Looking preference (s) \n <- Expected ------- Unexpected ->")+
  theme(axis.text.x = element_text(angle = 90, hjust=0.95, vjust =0.2), legend.position = "none")+
  scale_fill_brewer(palette = "Dark2")

plot2.constraints <- ggplot(ind.data %>% filter(task == "constraints"),
                 aes(condition, look_pref, fill=paper)) +
  geom_boxplot(outlier.alpha = 0.2, position = position_dodge2(preserve = "single"))+
  geom_hline(yintercept = 0)+
  geom_point(alpha=0.2)+
  stat_summary(geom="point", fun="mean", colour="black")+
  stat_summary(geom="errorbar", fun.data="mean_se", width=0.2, colour="black")+
  theme_cowplot(18)+
  facet_grid(~paper,scales = "free_x",space = "free", labeller = label_wrap_gen())+
  ylab("")+
  xlab("")+
  theme(axis.text.x = element_text(angle = 90, hjust=0.95, vjust =0.2), legend.position = "none")+
  scale_fill_manual(values = wes_palette("Royal2"))


plot2.goals

plot2.constraints

(plot2.goals | plot2.constraints) + plot_layout(widths = c(2, 1)) + plot_annotation(tag_levels = 'A')

(plot3 <- ggplot(ind.data,
                 aes(ageday, look_pref, colour=paper)) +
  # geom_boxplot(outlier.alpha = 0.2, position = position_dodge2(preserve = "single"))+
  geom_hline(yintercept = 0)+
  geom_point()+
  geom_smooth(method="lm")+
  # stat_summary(geom="point", fun="mean", colour="black")+
  # stat_summary(geom="errorbar", fun.data="mean_se", width=0.2, colour="black")+
  theme_cowplot(12)+
  facet_grid(~task+paper)+
  ylab("Looking preference (s) \n <- Expected ------- Unexpected ->")+
  theme(axis.text.x = element_text(angle = 90, hjust=0.95, vjust =0.2)))

# `geom_smooth()` using formula 'y ~ x'

Constraints Task

constraints <- ind.data %>% filter(task =="constraints")

constraints$actor_hand

#   [1] mittened mittened mittened mittened mittened mittened mittened mittened
#   [9] mittened mittened mittened mittened mittened mittened mittened mittened
#  [17] mittened mittened mittened mittened mittened mittened mittened mittened
#  [25] mittened mittened mittened mittened mittened mittened mittened mittened
#  [33] mittened mittened mittened mittened mittened mittened mittened mittened
#  [41] mittened mittened mittened mittened mittened mittened mittened mittened
#  [49] mittened mittened mittened mittened mittened mittened mittened mittened
#  [57] mittened mittened mittened mittened mittened mittened mittened mittened
#  [65] mittened mittened mittened mittened mittened mittened mittened mittened
#  [73] mittened mittened mittened mittened mittened mittened mittened mittened
#  [81] mittened mittened mittened mittened mittened mittened mittened mittened
#  [89] mittened mittened mittened mittened mittened mittened mittened mittened
#  [97] mittened mittened mittened mittened mittened mittened mittened mittened
# [105] mittened mittened mittened mittened mittened mittened mittened mittened
# [113] gloved   gloved   gloved   gloved   gloved   gloved   gloved   gloved  
# [121] gloved   gloved   gloved   gloved   gloved   gloved   gloved   gloved  
# [129] gloved   gloved   gloved   gloved   gloved   gloved   gloved   gloved  
# [137] gloved   gloved   gloved   gloved   gloved   gloved   gloved   gloved  
# [145] gloved   gloved   gloved   gloved   gloved   gloved   gloved   gloved  
# [153] gloved   gloved   gloved   gloved   gloved   gloved   gloved   gloved  
# [161] gloved   gloved   gloved   gloved   gloved   gloved   gloved   gloved  
# [169] gloved   gloved   gloved   gloved   gloved   gloved   gloved   gloved  
# [177] gloved   gloved   gloved   gloved   gloved   gloved   gloved   gloved  
# [185] gloved   gloved   gloved   gloved   gloved   gloved   gloved   gloved  
# [193] gloved   gloved   gloved   gloved   gloved   gloved   gloved   gloved  
# [201] gloved   gloved   gloved   gloved   gloved   gloved   gloved   gloved  
# [209] gloved   gloved   gloved   gloved   gloved   gloved   gloved   gloved  
# [217] gloved   gloved   gloved   gloved   gloved   gloved   gloved   gloved  
# [225] gloved   gloved   gloved   gloved   gloved   gloved   gloved   gloved  
# [233] gloved   gloved   gloved   gloved   gloved   gloved   gloved   gloved  
# [241] gloved   gloved   gloved   gloved   bare     bare     bare     bare    
# [249] bare     bare     bare     bare     bare     bare     bare     bare    
# [257] bare     bare     bare     bare     bare     bare     bare     bare    
# Levels: bare gloved mittened

constraints$action_consequence <- relevel(constraints$action_consequence, ref = "none")
constraints$agent_efficient_fam <- relevel(constraints$agent_efficient_fam, ref = "no")
constraints$training_yesno <- relevel(constraints$training_yesno, ref = "no")
constraints$action_causal <- relevel(constraints$action_causal, ref = "no")

Baseline Effect

constraints.baseline.data <- constraints %>%
  filter(condition %in% c("3_notraining",
                          "1_pickupglove",
                          "2_pickupbarehand"))
  
constraints.b1 <- lmer(data = constraints.baseline.data,
                       formula = look_pref ~ 1 + ageday + (1|condition))
summary(constraints.b1)

# Linear mixed model fit by REML. t-tests use Satterthwaite's method [
# lmerModLmerTest]
# Formula: look_pref ~ 1 + ageday + (1 | condition)
#    Data: constraints.baseline.data
# 
# REML criterion at convergence: 366.4
# 
# Scaled residuals: 
#     Min      1Q  Median      3Q     Max 
# -2.7189 -0.5297 -0.0435  0.7252  2.2686 
# 
# Random effects:
#  Groups    Name        Variance Std.Dev.
#  condition (Intercept)  2.03    1.42    
#  Residual              25.37    5.04    
# Number of obs: 60, groups:  condition, 3
# 
# Fixed effects:
#             Estimate Std. Error      df t value Pr(>|t|)
# (Intercept)  -2.3930     9.0407 57.2416   -0.26     0.79
# ageday        0.0258     0.0832 56.0050    0.31     0.76
# 
# Correlation of Fixed Effects:
#        (Intr)
# ageday -0.993

constraints.b1.table <- sjPlot:: tab_model(constraints.b1,
                                           show.std =TRUE,
                                           show.stat=TRUE,
                                           show.df=TRUE)

constraints.b1.beta <- summary(effectsize::standardize(constraints.b1))

constraints.baseline <- cbind(
  gen.beta(effectsize::standardize(constraints.b1)),
  gen.m(constraints.b1),
  gen.ci(constraints.b1)[3:4,]
)

# Computing profile confidence intervals ...

constraints.baseline

#                                beta        b      se    df       t      p
# (Intercept) -0.00000000000000002645 -2.39296 9.04071 57.24 -0.2647 0.7922
# ageday       0.03966247939082377660  0.02582 0.08316 56.01  0.3105 0.7574
#                lower   upper
# (Intercept) -20.1694 15.4724
# ageday       -0.1389  0.1897

cooks1 <- cooks.distance(constraints.b1, group = "subj")

# Warning in cooks.distance.lmerMod(constraints.b1, group = "subj"): group is not
# a valid argument for this function. As of version 0.4.0, group has been replaced
# by level.

dotplot_diag(x = cooks1, cutoff = "internal", name = "cooks.distance",  index=constraints.baseline.data$subj) + ylab("Cook's distance") + xlab("subjID")

excluded.subs <- c("90", "553", "86") 

constraints.b1.cooks <- lmer(data = constraints.baseline.data %>%
                         filter(!subj %in% excluded.subs),
                       formula = look_pref ~ 1 + ageday + (1|condition))

summary(constraints.b1.cooks)

# Linear mixed model fit by REML. t-tests use Satterthwaite's method [
# lmerModLmerTest]
# Formula: look_pref ~ 1 + ageday + (1 | condition)
#    Data: constraints.baseline.data %>% filter(!subj %in% excluded.subs)
# 
# REML criterion at convergence: 338.6
# 
# Scaled residuals: 
#     Min      1Q  Median      3Q     Max 
# -2.5544 -0.5380 -0.0717  0.7259  2.4971 
# 
# Random effects:
#  Groups    Name        Variance Std.Dev.
#  condition (Intercept)  1.14    1.07    
#  Residual              21.61    4.65    
# Number of obs: 57, groups:  condition, 3
# 
# Fixed effects:
#             Estimate Std. Error      df t value Pr(>|t|)
# (Intercept)  -9.5734     9.0731 54.2948   -1.06     0.30
# ageday        0.0926     0.0838 53.5855    1.11     0.27
# 
# Correlation of Fixed Effects:
#        (Intr)
# ageday -0.995

plot(allEffects(constraints.b1.cooks))

constraints.baseline.cooks <- cbind(
  gen.beta(effectsize::standardize(constraints.b1.cooks)),
  gen.m(constraints.b1.cooks),
  gen.ci(constraints.b1.cooks)[3:4,]
)

# Computing profile confidence intervals ...

For standard versions of the constraints task, we found no differential looking between the expected and unexpected events, Mean_diff = 0.395 seconds, [-0.139,0.19], ß=0.04, B=0.026, SE=0.083, p=0.757, two-tailed. This result held when we removed 3 influential observations, identified using Cook’s distance, [-0.07,0.263], ß=0.145, B=0.093, SE=0.084, p=0.274, two-tailed.

Intervention Analysis

constraints.m1 <- lmer(data = constraints,
     formula = look_pref ~ training_yesno + action_causal + action_consequence + actor_hand + agent_efficient_fam + ageday + (1|condition) + (1|experiment) + (1|paper))

# boundary (singular) fit: see ?isSingular

# Warning in as_lmerModLT(model, devfun): Model may not have converged with 1
# eigenvalue close to zero: -3.1e-11

summary(constraints.m1)

# Linear mixed model fit by REML. t-tests use Satterthwaite's method [
# lmerModLmerTest]
# Formula: look_pref ~ training_yesno + action_causal + action_consequence +  
#     actor_hand + agent_efficient_fam + ageday + (1 | condition) +  
#     (1 | experiment) + (1 | paper)
#    Data: constraints
# 
# REML criterion at convergence: 1682
# 
# Scaled residuals: 
#    Min     1Q Median     3Q    Max 
# -3.729 -0.510 -0.059  0.561  3.866 
# 
# Random effects:
#  Groups     Name        Variance Std.Dev.
#  condition  (Intercept)  0.000   0.000   
#  experiment (Intercept)  0.000   0.000   
#  paper      (Intercept)  0.121   0.347   
#  Residual               35.324   5.943   
# Number of obs: 264, groups:  condition, 12; experiment, 5; paper, 2
# 
# Fixed effects:
#                      Estimate Std. Error       df t value Pr(>|t|)    
# (Intercept)           -2.8929     4.6171 256.0000   -0.63  0.53151    
# training_yesno1       -2.1859     0.6174 256.0000   -3.54  0.00047 ***
# action_causal1        -2.3185     0.5944 256.0000   -3.90  0.00012 ***
# action_consequence1   -1.0693     0.7763 256.0000   -1.38  0.16959    
# actor_hand1            1.3018     1.0518 256.0000    1.24  0.21698    
# actor_hand2            0.8084     1.0518 256.0000    0.77  0.44282    
# agent_efficient_fam1  -1.7116     0.5377 256.0000   -3.18  0.00164 ** 
# ageday                 0.0231     0.0419 256.0000    0.55  0.58139    
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# Correlation of Fixed Effects:
#             (Intr) trnn_1 actn_cs1 actn_cn1 actr_1 actr_2 agn__1
# tranng_ysn1 -0.062                                              
# action_csl1  0.143  0.085                                       
# actn_cnsqn1  0.010  0.065  0.349                                
# actor_hand1  0.124 -0.227  0.001   -0.360                       
# actor_hand2 -0.060 -0.227  0.000    0.721   -0.597              
# agnt_ffcn_1  0.092  0.314  0.275    0.210    0.000  0.000       
# ageday      -0.982  0.017 -0.053   -0.036   -0.010 -0.008 -0.018
# convergence code: 0
# boundary (singular) fit: see ?isSingular

resid_panel(constraints.m1, plots = "all",
                          smoother = TRUE,
                          qqbands = TRUE)

# `geom_smooth()` using formula 'y ~ x'

# `geom_smooth()` using formula 'y ~ x'

cooks1 <- cooks.distance(constraints.m1, group = "subj")

# Warning in cooks.distance.lmerMod(constraints.m1, group = "subj"): group is not
# a valid argument for this function. As of version 0.4.0, group has been replaced
# by level.

dotplot_diag(x = cooks1, cutoff = "internal", name = "cooks.distance", index=constraints$subj) + ylab("Cook's distance") + xlab("ID")

sjPlot::plot_model(constraints.m1,
                   type = "est",
                   colors = "bw",
                   sort.est=TRUE,
                   axis.title=c("Effect (Unexpected - Expected in seconds)"),
                   show.values=TRUE,
                   show.p=TRUE)

constraints.m1.table <- sjPlot:: tab_model(constraints.m1,
                                           show.std =TRUE,
                                           show.stat=TRUE,
                                           show.df=TRUE)

# boundary (singular) fit: see ?isSingular

constraints.m1.beta <- summary(effectsize::standardize(constraints.m1))

# boundary (singular) fit: see ?isSingular

constraints.interventions <- cbind(
  gen.beta(effectsize::standardize(constraints.m1)),
  gen.m(constraints.m1),
  gen.ci(constraints.m1)[3:10,]
)

# boundary (singular) fit: see ?isSingular

# Computing profile confidence intervals ...

constraints.cooks<- constraints[which(cooks1 <= 4/264),]

# where are the influential observations?
nrow(constraints)

# [1] 264

nrow(constraints.cooks)

# [1] 249

constraints.summary <- constraints %>%
  group_by(paper,condition) %>%
  summarise(n = n())

# `summarise()` regrouping output by 'paper' (override with `.groups` argument)

constraints.cooksn <- constraints.cooks %>%
  group_by(paper,condition) %>%
  summarise(n_cooks = n())

# `summarise()` regrouping output by 'paper' (override with `.groups` argument)

constraints.where.influential <- full_join(constraints.summary,constraints.cooksn)

# Joining, by = c("paper", "condition")

constraints.m1.cooks <- lmer(data = constraints.cooks,
     formula = look_pref ~ training_yesno + action_causal + action_consequence + actor_hand + agent_efficient_fam + ageday + (1|condition) + (1|experiment) + (1|paper),
     REML=FALSE)

# boundary (singular) fit: see ?isSingular

summary(constraints.m1.cooks)

# Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
#   method [lmerModLmerTest]
# Formula: look_pref ~ training_yesno + action_causal + action_consequence +  
#     actor_hand + agent_efficient_fam + ageday + (1 | condition) +  
#     (1 | experiment) + (1 | paper)
#    Data: constraints.cooks
# 
#      AIC      BIC   logLik deviance df.resid 
#   1484.8   1527.0   -730.4   1460.8      237 
# 
# Scaled residuals: 
#     Min      1Q  Median      3Q     Max 
# -2.9748 -0.6479 -0.0429  0.6569  3.1180 
# 
# Random effects:
#  Groups     Name        Variance Std.Dev.
#  condition  (Intercept)  0.0     0.00    
#  experiment (Intercept)  0.0     0.00    
#  paper      (Intercept)  0.0     0.00    
#  Residual               20.7     4.55    
# Number of obs: 249, groups:  condition, 12; experiment, 5; paper, 2
# 
# Fixed effects:
#                      Estimate Std. Error       df t value   Pr(>|t|)    
# (Intercept)           -4.0078     3.6877 249.0000   -1.09      0.278    
# training_yesno1       -1.9535     0.4773 249.0000   -4.09 0.00005773 ***
# action_causal1        -2.5634     0.4768 249.0000   -5.38 0.00000017 ***
# action_consequence1   -1.2255     0.6202 249.0000   -1.98      0.049 *  
# actor_hand1            0.7581     0.8187 249.0000    0.93      0.355    
# actor_hand2            0.9867     0.8311 249.0000    1.19      0.236    
# agent_efficient_fam1  -1.9383     0.4261 249.0000   -4.55 0.00000844 ***
# ageday                 0.0281     0.0333 249.0000    0.84      0.400    
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# Correlation of Fixed Effects:
#             (Intr) trnn_1 actn_cs1 actn_cn1 actr_1 actr_2 agn__1
# tranng_ysn1 -0.061                                              
# action_csl1  0.161  0.076                                       
# actn_cnsqn1  0.007  0.058  0.330                                
# actor_hand1  0.120 -0.226  0.001   -0.377                       
# actor_hand2 -0.047 -0.223  0.002    0.744   -0.644              
# agnt_ffcn_1  0.124  0.313  0.248    0.190    0.000  0.001       
# ageday      -0.983  0.018 -0.070   -0.035   -0.009 -0.025 -0.050
# convergence code: 0
# boundary (singular) fit: see ?isSingular

constraints.m1.cooks.beta <- lmer(data = constraints.cooks,
     formula = scale(look_pref) ~ training_yesno + action_causal + action_consequence + actor_hand + agent_efficient_fam + scale(ageday) + (1|condition) + (1|experiment) + (1|paper),
     REML=FALSE)

# boundary (singular) fit: see ?isSingular

summary(constraints.m1.cooks.beta)

# Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
#   method [lmerModLmerTest]
# Formula: 
# scale(look_pref) ~ training_yesno + action_causal + action_consequence +  
#     actor_hand + agent_efficient_fam + scale(ageday) + (1 | condition) +  
#     (1 | experiment) + (1 | paper)
#    Data: constraints.cooks
# 
#      AIC      BIC   logLik deviance df.resid 
#    683.5    725.7   -329.8    659.5      237 
# 
# Scaled residuals: 
#     Min      1Q  Median      3Q     Max 
# -2.9748 -0.6479 -0.0429  0.6569  3.1180 
# 
# Random effects:
#  Groups     Name        Variance Std.Dev.
#  condition  (Intercept) 0.000    0.00    
#  experiment (Intercept) 0.000    0.00    
#  paper      (Intercept) 0.000    0.00    
#  Residual               0.828    0.91    
# Number of obs: 249, groups:  condition, 12; experiment, 5; paper, 2
# 
# Fixed effects:
#                      Estimate Std. Error       df t value   Pr(>|t|)    
# (Intercept)           -0.3453     0.1338 249.0000   -2.58      0.010 *  
# training_yesno1       -0.3908     0.0955 249.0000   -4.09 0.00005773 ***
# action_causal1        -0.5129     0.0954 249.0000   -5.38 0.00000017 ***
# action_consequence1   -0.2452     0.1241 249.0000   -1.98      0.049 *  
# actor_hand1            0.1517     0.1638 249.0000    0.93      0.355    
# actor_hand2            0.1974     0.1663 249.0000    1.19      0.236    
# agent_efficient_fam1  -0.3878     0.0853 249.0000   -4.55 0.00000844 ***
# scale(ageday)          0.0490     0.0580 249.0000    0.84      0.400    
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# Correlation of Fixed Effects:
#             (Intr) trnn_1 actn_cs1 actn_cn1 actr_1 actr_2 agn__1
# tranng_ysn1 -0.238                                              
# action_csl1  0.510  0.076                                       
# actn_cnsqn1 -0.151  0.058  0.330                                
# actor_hand1  0.613 -0.226  0.001   -0.377                       
# actor_hand2 -0.393 -0.223  0.002    0.744   -0.644              
# agnt_ffcn_1  0.415  0.313  0.248    0.190    0.000  0.001       
# scale(agdy) -0.058  0.018 -0.070   -0.035   -0.009 -0.025 -0.050
# convergence code: 0
# boundary (singular) fit: see ?isSingular

plot_model(constraints.m1.cooks.beta)

sjPlot::plot_model(constraints.m1.cooks,
                   type = "pred",
                   # colors = "bw",
                   sort.est=TRUE,
                   axis.title=c("Effect (Unexpected - Expected in seconds)"),
                   show.values=TRUE,
                   show.p=TRUE)

# $training_yesno

# 
# $action_causal

# 
# $action_consequence

# 
# $actor_hand

# 
# $agent_efficient_fam

# 
# $ageday

plot(allEffects(constraints.m1.cooks))

sjPlot::plot_model(constraints.m1.cooks,
                   type = "est",
                   colors = "bw",
                   sort.est=TRUE,
                    axis.labels= rev(c("mittened", "gloved", "age", "state change", "actor initially rational", "sticky mittens", "action effect on contact")),
                   axis.title=c("Effect (Unexpected - Expected in seconds)"),
                   show.values=TRUE,
                   show.p=TRUE)

constraints.interventions.cooks <- cbind(
  gen.beta(effectsize::standardize(constraints.m1.cooks)),
  gen.m(constraints.m1.cooks),
  gen.ci(constraints.m1.cooks)[3:10,]
)

# boundary (singular) fit: see ?isSingular
# Computing profile confidence intervals ...

When we compared the effects of different variants we found that motor training increased infants’ looking preference for the unexpected event ([4.175,4.978], ß=-0.391, B=-1.953, SE=0.477, p<.001, two-tailed). We also found that other interventions on the actions infants see, including seeing an action that causes an effect on contact ([-11.263,3.247], ß=-0.513, B=-2.563, SE=0.477, p<.001, two-tailed), and seeing a reaching action that results in an object changing state, rather than picking up that object ([-2.892,-1.014], ß=-0.245, B=-1.226, SE=0.62, p=0.049, two-tailed). These analyses included 249/264 total infants - the other 15 were classified as influential using Cook’s Distance. Including all subjects in the analysis produced similar results, except that the state change effect crossed the threshold into non-significance ([-3.382,-0.989], ß=-0.173, B=-1.069, SE=0.776, p=0.17, two-tailed). We did not find an effect of age, [-0.649,2.623], ß=0.049, B=0.028, SE=0.033, p=0.4, two-tailed.

Goals Task

Baseline Effect

goals <- ind.data %>% filter(task =="goals")

goals.baseline.data <- goals %>%
  filter(condition %in% c("1_control",
                          "1_twoobject") &
           paper %in% c("gerson2014a",
                        "luo2011"))
  
goals.b1 <- lmer(data = goals.baseline.data,
                       formula = look_pref ~ 1 + ageday + (1|condition))

goals.b1.table <- sjPlot:: tab_model(goals.b1,
                                           show.std =TRUE,
                                           show.stat=TRUE,
                                           show.df=TRUE)

goals.b1.beta <- summary(effectsize::standardize(goals.b1))

goals.baseline <- cbind(
  gen.beta(effectsize::standardize(goals.b1)),
  gen.m(goals.b1),
  gen.ci(goals.b1)[3:4,]
)

# Computing profile confidence intervals ...

cooks1 <- cooks.distance(goals.b1, group = "subj")

# Warning in cooks.distance.lmerMod(goals.b1, group = "subj"): group is not a
# valid argument for this function. As of version 0.4.0, group has been replaced
# by level.

dotplot_diag(x = cooks1, cutoff = "internal", name = "cooks.distance",  index=goals.baseline.data$subj) + ylab("Cook's distance") + xlab("subjID")

excluded.subs <- c("2", "7", "6", "11") 

goals.b1.cooks <- lmer(data = goals.baseline.data %>%
                         filter(!subj %in% excluded.subs),
                       formula = look_pref ~ 1 + ageday + (1|condition))

summary(goals.b1.cooks)

# Linear mixed model fit by REML. t-tests use Satterthwaite's method [
# lmerModLmerTest]
# Formula: look_pref ~ 1 + ageday + (1 | condition)
#    Data: goals.baseline.data %>% filter(!subj %in% excluded.subs)
# 
# REML criterion at convergence: 239.2
# 
# Scaled residuals: 
#     Min      1Q  Median      3Q     Max 
# -1.6626 -0.5931  0.0995  0.4841  2.8202 
# 
# Random effects:
#  Groups    Name        Variance Std.Dev.
#  condition (Intercept) 202      14.2    
#  Residual              106      10.3    
# Number of obs: 32, groups:  condition, 2
# 
# Fixed effects:
#             Estimate Std. Error     df t value Pr(>|t|)
# (Intercept)   44.720     27.063 21.547    1.65     0.11
# ageday        -0.355      0.236 29.073   -1.51     0.14
# 
# Correlation of Fixed Effects:
#        (Intr)
# ageday -0.925

goals.baseline.cooks <- cbind(
  gen.beta(effectsize::standardize(goals.b1.cooks)),
  gen.m(goals.b1.cooks),
  gen.ci(goals.b1.cooks)[3:4,]
)

# Computing profile confidence intervals ...

# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep

# Warning in profile.merMod(object, which = parm, signames = oldNames, ...): non-
# monotonic profile for (Intercept)

# Warning in confint.thpr(pp, level = level, zeta = zeta): bad spline fit for
# (Intercept): falling back to linear interpolation

For standard versions of the goals task, we found differential looking between the expected and unexpected events, Mean_diff = 3.918 seconds, [-0.931,-0.081], ß=-0.318, B=-0.503, SE=0.214, p=0.025, two-tailed. However, this was driven by 4 influential observations - without these 4 observations, there was no longer a reliable looking preference, [-0.815,0.126], ß=-0.21, B=-0.355, SE=0.236, p=0.143, two-tailed.

Intervention Analysis

goals.m1 <- lmer(data = goals,
     formula = look_pref ~ training_yesno  + action_consequence + location_object_goal_ambiguous + agent + bothobjects_present_visible_fam + ageday + (1|condition) + (1|experiment) + (1|paper),
     control = lmerControl(optimizer="Nelder_Mead"))

# boundary (singular) fit: see ?isSingular

summary(goals.m1)

# Linear mixed model fit by REML. t-tests use Satterthwaite's method [
# lmerModLmerTest]
# Formula: 
# look_pref ~ training_yesno + action_consequence + location_object_goal_ambiguous +  
#     agent + bothobjects_present_visible_fam + ageday + (1 | condition) +  
#     (1 | experiment) + (1 | paper)
#    Data: goals
# Control: lmerControl(optimizer = "Nelder_Mead")
# 
# REML criterion at convergence: 2842
# 
# Scaled residuals: 
#    Min     1Q Median     3Q    Max 
# -4.049 -0.458 -0.028  0.416  2.921 
# 
# Random effects:
#  Groups     Name        Variance Std.Dev.
#  condition  (Intercept)  19.6     4.43   
#  paper      (Intercept)   0.0     0.00   
#  experiment (Intercept)   0.0     0.00   
#  Residual               153.4    12.39   
# Number of obs: 362, groups:  condition, 13; paper, 6; experiment, 2
# 
# Fixed effects:
#                                  Estimate Std. Error       df t value Pr(>|t|)
# (Intercept)                       10.3237     6.2887  96.2681    1.64    0.104
# training_yesno1                   -1.9519     2.4419   4.8329   -0.80    0.462
# action_consequence1                2.8811     2.4462   7.9680    1.18    0.273
# location_object_goal_ambiguous1    5.1289     2.4944   8.5843    2.06    0.071
# agent1                            -5.4390     2.9960   6.3026   -1.82    0.117
# agent2                             6.4666     2.7646  19.3741    2.34    0.030
# bothobjects_present_visible_fam1  -5.3013     2.0586  15.0765   -2.58    0.021
# ageday                            -0.0733     0.0531 352.2521   -1.38    0.168
#                                   
# (Intercept)                       
# training_yesno1                   
# action_consequence1               
# location_object_goal_ambiguous1  .
# agent1                            
# agent2                           *
# bothobjects_present_visible_fam1 *
# ageday                            
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# Correlation of Fixed Effects:
#             (Intr) trnn_1 actn_1 lc___1 agent1 agent2 bt___1
# tranng_ysn1 -0.295                                          
# actn_cnsqn1 -0.246  0.002                                   
# lctn_bjc__1  0.067  0.002  0.482                            
# agent1       0.011  0.547 -0.169  0.111                     
# agent2      -0.013 -0.297  0.024 -0.170 -0.796              
# bthbjct___1  0.119  0.002 -0.245 -0.170  0.387 -0.563       
# ageday      -0.875  0.043  0.053  0.038 -0.042  0.017  0.057
# convergence code: 0
# boundary (singular) fit: see ?isSingular

goals %>%
  select(paper, experiment, condition, training_yesno, object_diff_size_huge,action_consequence,location_object_goal_ambiguous,agent,bothobjects_present_visible_fam) %>%
  distinct()

#          paper experiment                  condition training_yesno
# 1      luo2011          1                1_twoobject             no
# 2      luo2011          1                1_oneobject             no
# 3      luo2011          2       2_differentpositions             no
# 4  gerson2014a          1                   1_active            yes
# 5  gerson2014a          1              1_observation             no
# 6  gerson2014a          1                  1_control             no
# 7  gerson2014b          1                   1_active            yes
# 8  gerson2014b          1              1_observation             no
# 9  gerson2014b          2           2_generalization            yes
# 10  choi_unpub          1              1_equidistant             no
# 11  choi_unpub          1               1_neartarget             no
# 12  choi_unpub          1                1_fartarget             no
# 13    choi2018          1                1_twoobject             no
# 14    choi2018          1                1_oneobject             no
# 15    choi2018          1                   1_hidden             no
# 16   woo_unpub          1   1_statechange_objectswap             no
# 17   woo_unpub          2 2_disambiguatingobjectgoal             no
#    object_diff_size_huge action_consequence location_object_goal_ambiguous
# 1                     no               none                            yes
# 2                     no               none                            yes
# 3                     no               none                             no
# 4                     no               none                            yes
# 5                     no               none                            yes
# 6                     no               none                            yes
# 7                     no               none                            yes
# 8                     no               none                            yes
# 9                     no               none                            yes
# 10                   yes               none                            yes
# 11                   yes               none                            yes
# 12                   yes               none                            yes
# 13                   yes               none                            yes
# 14                   yes               none                            yes
# 15                   yes               none                            yes
# 16                    no       state_change                            yes
# 17                    no       state_change                             no
#     agent bothobjects_present_visible_fam
# 1  object                             yes
# 2  object                              no
# 3  object                              no
# 4    hand                             yes
# 5    hand                             yes
# 6    hand                             yes
# 7    hand                             yes
# 8    hand                             yes
# 9    hand                             yes
# 10 person                             yes
# 11 person                             yes
# 12 person                              no
# 13 person                             yes
# 14 person                             yes
# 15 person                             yes
# 16 person                             yes
# 17 person                             yes

# originally pre-registered model was rank deficient, needed to drop 2 variables
# dropped action_causal because it is redundant with action_consequence (i.e. there are no non-causal actions that are state changes)
# dropped object_diff_size_huge also because it is almost completely redundant with agent (i.e. all studies with a full human agent also are "yes" for "object_diff_size_huge")

# aa <- allFit(goals.m1)

summary(goals.m1)

# Linear mixed model fit by REML. t-tests use Satterthwaite's method [
# lmerModLmerTest]
# Formula: 
# look_pref ~ training_yesno + action_consequence + location_object_goal_ambiguous +  
#     agent + bothobjects_present_visible_fam + ageday + (1 | condition) +  
#     (1 | experiment) + (1 | paper)
#    Data: goals
# Control: lmerControl(optimizer = "Nelder_Mead")
# 
# REML criterion at convergence: 2842
# 
# Scaled residuals: 
#    Min     1Q Median     3Q    Max 
# -4.049 -0.458 -0.028  0.416  2.921 
# 
# Random effects:
#  Groups     Name        Variance Std.Dev.
#  condition  (Intercept)  19.6     4.43   
#  paper      (Intercept)   0.0     0.00   
#  experiment (Intercept)   0.0     0.00   
#  Residual               153.4    12.39   
# Number of obs: 362, groups:  condition, 13; paper, 6; experiment, 2
# 
# Fixed effects:
#                                  Estimate Std. Error       df t value Pr(>|t|)
# (Intercept)                       10.3237     6.2887  96.2681    1.64    0.104
# training_yesno1                   -1.9519     2.4419   4.8329   -0.80    0.462
# action_consequence1                2.8811     2.4462   7.9680    1.18    0.273
# location_object_goal_ambiguous1    5.1289     2.4944   8.5843    2.06    0.071
# agent1                            -5.4390     2.9960   6.3026   -1.82    0.117
# agent2                             6.4666     2.7646  19.3741    2.34    0.030
# bothobjects_present_visible_fam1  -5.3013     2.0586  15.0765   -2.58    0.021
# ageday                            -0.0733     0.0531 352.2521   -1.38    0.168
#                                   
# (Intercept)                       
# training_yesno1                   
# action_consequence1               
# location_object_goal_ambiguous1  .
# agent1                            
# agent2                           *
# bothobjects_present_visible_fam1 *
# ageday                            
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# Correlation of Fixed Effects:
#             (Intr) trnn_1 actn_1 lc___1 agent1 agent2 bt___1
# tranng_ysn1 -0.295                                          
# actn_cnsqn1 -0.246  0.002                                   
# lctn_bjc__1  0.067  0.002  0.482                            
# agent1       0.011  0.547 -0.169  0.111                     
# agent2      -0.013 -0.297  0.024 -0.170 -0.796              
# bthbjct___1  0.119  0.002 -0.245 -0.170  0.387 -0.563       
# ageday      -0.875  0.043  0.053  0.038 -0.042  0.017  0.057
# convergence code: 0
# boundary (singular) fit: see ?isSingular

sjPlot::plot_model(goals.m1,
                   type = "est",
                   colors = "bw",
                   sort.est=TRUE,
                   axis.title=c("Effect (Unexpected - Expected in seconds)"),
                   show.values=TRUE,
                   show.p=TRUE)

sjPlot:: tab_model(goals.m1, show.stat=TRUE)

	look pref
Predictors	Estimates	CI	Statistic	p
(Intercept)	10.32	-2.00 – 22.65	1.64	0.101
training_yesno1	-1.95	-6.74 – 2.83	-0.80	0.424
action_consequence1	2.88	-1.91 – 7.68	1.18	0.239
location_object_goal_ambiguous1	5.13	0.24 – 10.02	2.06	0.040
agent1	-5.44	-11.31 – 0.43	-1.82	0.069
agent2	6.47	1.05 – 11.89	2.34	0.019
bothobjects_present_visible_fam1	-5.30	-9.34 – -1.27	-2.58	0.010
ageday	-0.07	-0.18 – 0.03	-1.38	0.168
Random Effects
σ²	153.43
τ₀₀ _condition	19.60
τ₀₀ _paper	0.00
τ₀₀ _experiment	0.00
N _condition	13
N _experiment	2
N _paper	6
Observations	362
Marginal R² / Conditional R²	0.104 / NA

plot(allEffects(goals.m1))

goals.interventions<- cbind(
  gen.beta(effectsize::standardize(goals.m1)),
  gen.m(goals.m1),
  gen.ci(goals.m1)[4:11,]
)

# boundary (singular) fit: see ?isSingular

# Computing profile confidence intervals ...

# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep

# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep

# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep

# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep

# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep

# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep

# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep

# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep

# Warning in profile.merMod(object, which = parm, signames = oldNames, ...): non-
# monotonic profile for (Intercept)

# Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
# profile: using minstep

# Warning in profile.merMod(object, which = parm, signames = oldNames, ...): non-
# monotonic profile for agent2

# Warning in confint.thpr(pp, level = level, zeta = zeta): bad spline fit for
# (Intercept): falling back to linear interpolation

# Warning in confint.thpr(pp, level = level, zeta = zeta): bad spline fit for
# agent2: falling back to linear interpolation

cooks1 <- cooks.distance(goals.m1, group = "subj")

# Warning in cooks.distance.lmerMod(goals.m1, group = "subj"): group is not a
# valid argument for this function. As of version 0.4.0, group has been replaced
# by level.

dotplot_diag(x = cooks1, cutoff = "internal", name = "cooks.distance",  index=goals$subj) + ylab("Cook's distance") + xlab("subjID")

goals.cooks <- goals[which(cooks1 <= 4/362),]


# where are the influential observations?
nrow(goals)

# [1] 362

nrow(goals.cooks)

# [1] 345

goals.summary <- goals %>%
  group_by(paper,condition) %>%
  summarise(n = n())

# `summarise()` regrouping output by 'paper' (override with `.groups` argument)

goals.cooksn <- goals.cooks %>%
  group_by(paper,condition) %>%
  summarise(n_cooks = n())

# `summarise()` regrouping output by 'paper' (override with `.groups` argument)

goals.where.influential <- full_join(goals.summary, goals.cooksn)

# Joining, by = c("paper", "condition")

goals.m1.cooks <- lmer(data = goals.cooks,
     formula = look_pref ~ training_yesno  + action_consequence + location_object_goal_ambiguous + agent + bothobjects_present_visible_fam + ageday + (1|condition) + (1|experiment) + (1|paper),
     control = lmerControl(optimizer="Nelder_Mead"))

# boundary (singular) fit: see ?isSingular

fixef(goals.m1.cooks)

#                      (Intercept)                  training_yesno1 
#                          6.06716                         -1.83931 
#              action_consequence1  location_object_goal_ambiguous1 
#                          2.33221                          5.15730 
#                           agent1                           agent2 
#                         -4.46125                          5.60015 
# bothobjects_present_visible_fam1                           ageday 
#                         -3.07025                         -0.01669

fixef(goals.m1.cooks,add.dropped=TRUE)

#                      (Intercept)                  training_yesno1 
#                          6.06716                         -1.83931 
#              action_consequence1  location_object_goal_ambiguous1 
#                          2.33221                          5.15730 
#                           agent1                           agent2 
#                         -4.46125                          5.60015 
# bothobjects_present_visible_fam1                           ageday 
#                         -3.07025                         -0.01669

summary(goals.m1.cooks)

# Linear mixed model fit by REML. t-tests use Satterthwaite's method [
# lmerModLmerTest]
# Formula: 
# look_pref ~ training_yesno + action_consequence + location_object_goal_ambiguous +  
#     agent + bothobjects_present_visible_fam + ageday + (1 | condition) +  
#     (1 | experiment) + (1 | paper)
#    Data: goals.cooks
# Control: lmerControl(optimizer = "Nelder_Mead")
# 
# REML criterion at convergence: 2595
# 
# Scaled residuals: 
#    Min     1Q Median     3Q    Max 
# -3.484 -0.459 -0.042  0.468  3.371 
# 
# Random effects:
#  Groups     Name        Variance Std.Dev.
#  condition  (Intercept)  13.3     3.64   
#  paper      (Intercept)   0.0     0.00   
#  experiment (Intercept)   0.0     0.00   
#  Residual               110.0    10.49   
# Number of obs: 345, groups:  condition, 13; paper, 6; experiment, 2
# 
# Fixed effects:
#                                  Estimate Std. Error       df t value Pr(>|t|)
# (Intercept)                        6.0672     5.4877 109.6470    1.11    0.271
# training_yesno1                   -1.8393     2.0190   4.7758   -0.91    0.406
# action_consequence1                2.3322     2.0883   8.2716    1.12    0.295
# location_object_goal_ambiguous1    5.1573     2.2202   9.8410    2.32    0.043
# agent1                            -4.4613     2.5736   6.7760   -1.73    0.128
# agent2                             5.6002     2.6245  23.8723    2.13    0.043
# bothobjects_present_visible_fam1  -3.0702     1.8673  12.4200   -1.64    0.125
# ageday                            -0.0167     0.0465 336.0909   -0.36    0.720
#                                   
# (Intercept)                       
# training_yesno1                   
# action_consequence1               
# location_object_goal_ambiguous1  *
# agent1                            
# agent2                           *
# bothobjects_present_visible_fam1  
# ageday                            
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# Correlation of Fixed Effects:
#             (Intr) trnn_1 actn_1 lc___1 agent1 agent2 bt___1
# tranng_ysn1 -0.284                                          
# actn_cnsqn1 -0.226  0.002                                   
# lctn_bjc__1  0.067  0.002  0.519                            
# agent1      -0.005  0.527 -0.154  0.123                     
# agent2       0.027 -0.260 -0.007 -0.201 -0.796              
# bthbjct___1  0.114  0.003 -0.247 -0.151  0.411 -0.545       
# ageday      -0.880  0.046  0.048  0.034 -0.023 -0.014  0.073
# convergence code: 0
# boundary (singular) fit: see ?isSingular

plot(allEffects(goals.m1.cooks))

goals.interventions.cooks<- cbind(
  gen.beta(effectsize::standardize(goals.m1.cooks)),
  gen.m(goals.m1.cooks),
  gen.ci(goals.m1.cooks)[4:11,]
)

# boundary (singular) fit: see ?isSingular
# Computing profile confidence intervals ...

We note that the analysis we reported deviates from our pre-registration plan. In that plan, we included 8 fixed effects, including a fixed effect for age in days, one picking out a control condition, and 6 others. However, 2 of these fixed effects were either completely or partially redundant with other predictors in the model, which caused a rank deficiency issue. Thus, we dropped these two fixed effects from the model. We report this new model below.

When we compared the effects of different variants, we found that making the goal of the agent’s actions unambiguous ([-1.191,5.201], ß=0.465, B=5.157, SE=2.22, p=0.043, two-tailed), or presenting a self-propelled object as an agent ([-8.285,-0.729], ß=0.505, B=5.6, SE=2.625, p=0.043, two-tailed), relative to a hand or a person, increased infants’ looking preference for the unexpected event. We did not find any other effects, including those of motor experience, or action consequences, when taking into account all of these predictors in the same model. These analyses included 345/362 total infants - the other 17 were classified as influential using Cook’s Distance. Including all subjects in the analysis produced similar results, except that the goal ambiguity effect crossed the threshold into marginal-significance ([-3.382,-0.989], ß=-0.173, B=-1.069, SE=0.776, p=0.17, two-tailed). We did not find an effect of age, [-5.763,0.015], ß=-0.019, B=-0.017, SE=0.047, p=0.72, two-tailed.

Meta-analysis 1

Shari Liu

3/24/2021

Summary of data sources

Plots

Constraints Task

Baseline Effect

Intervention Analysis

Goals Task

Baseline Effect

Intervention Analysis